AD-A208  725 


Systems 

Optimization 

Laboratory 


An  Efficient  Exact  Algorithm  for  the 
’’LEAST  SQUARES”  Image  Registration  Problem  * 

by 

Karel  Zikan 

TECHNICAL  REPORT  SOL 
May  1989 


DTiC 


ELECTE 
JUNO  6  19891  i 


Department  of  Operations  Research 
Stanford  University 
Stanford,  CA  94305 

“DISTmuflO'N  STAfEMOrni  I 


Appro7«d  tax  publle  mIaom; 
Dtatrlbution  Umlimlted 


SYSTEMS  OPTIMIZATION  LABORATORY 
DEPARTMENT  OF  OPERATIONS  RESEARCH 
STANFORD  UNIVERSITY 
STANFORD,  CALIFORNIA  94305-4022 


An  EfRciciit  Exact  Algorithm  for  the 
’’LEAST  SQUARES”  Image  Registration  Problem  * 

by 

Karel  Zikan 

TECHNICAL  REPORT  SOL  89-5 
May  1989 


DTIC 

ELECTE 
JUNO  6  1989 

H 


0 


Reoearch  and  reproduction  of  this  report  were  partially  supported  by  the  National  Science  Founda¬ 
tion  Grant  DMS-8800589;  U  S.  Departn:ent  of  Energy  Grant  DE-’"G03-87ER25028,  and  Office  of 
Naval  Research  Grant  N00014-89-J- 1659. 

Any  oiMiHons.  findings,  and  conclusions  or  recommendations  expressed  in  thi.=  |uiblicalion  are  those 
of  the  author(s)  and  do  NOT  necessarily  reflect  the  views  of  the  above  sponsors. 

Reproduction  m  whole  or  in  part  is  permitted  for  any  purposes  of  the  United  States  Government. 
Ihis  document  has  been  approved  for  public  release  and  sale,  its  distribution  is  unlimited. 


This  paper  won  Honorable  Mention  in  the  1989  Nicholson  Student  Paper  Competition,  sponsored 
by  the  Operations  Research  Society  of  America. 


AN  EFFICIENT  EXACT  ALGORITHM  FOR  THE 
LEAST-SQUARES”  IMAGE  REGISTRATION  PROBLEM 


Karel  ZIKAN 


Department  of  Operations  Research 
STANFORD  UNIVERSITY 
Stanford,  CA  94305 

and 


Hughes  Artificial  Intelligence  Center 
Hughes  Research  Laboratories 
3011  Malibu  Cwyon  Rd.,  Malibu,  Ca  90265 


Abstract 

Image  registration  involves  estimating  how  one  set  of  n-dimensional  points  is  rotated,  scaled,  and 
translated  into  a  second  set  of  n-dimensional  points.  In  practice,  n  is  usually  2  or  3.  We  give  an  exact 
algorithm  to  solve  the  "least-squares’*  formulation  of  the  two-dimensional  registration  problem.  The 
algorithm,  which  is  based  on  parametric  linear  programming,  can  be  viewed  as  a  refinement  of  the 
0{k^)  approximation  method  proposed  by  Zikan  and  Silberberg^(l3].  The  approach  can  be  extended 
to  handle  registration  of  images  of  different  cardinalities. 
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1.  Background  of  the  Problem 


The  basic  image  registration  problem  is  stated  in  zikaN  and  silberberg  [13];  “Assume  thet  a 
collection  of  n-dimensional  points  undergoes  an  affine  transformation  made  up  of  a  rotation,  a  scale 
change,  and  a  translation;  moreover,  assume  that  noise  (random  and/or  systematic)  is  added  to 
each  transformed  point.  Knowing  the  positions  of  the  original  and  transformed  points,  but  not  their 
identities,  can  the  trsms formation,  the  noise,  and  the  point-to-point  matching  be  recovered?" 

The  image  registration  problem  is  fundamental  in  robotics,  as  the  ability  to  “register”  images  is  a 
necessary  prerequisite  for  “on-board”  automated  visual  reasoning  of  mobile  agents.  Many  research 
papers  liave  been  written  on  the  subject.  Henry  S.  Baird  [4]  (ACM  Distinguished  Thesis  Series, 
19S4)  outlines  some  of  the  older,  mostly  heuristic  approaches.  Many  of  the  methods  employ  the 
so-called  pri  ned  tree  approach;  see  e.g.,  gennery  [6]  and  wong  &nd  salay  [12],  where  the  tree  of 
partial  matchings  is  searched  for  the  desired  solution.  All  permutations  of  the  image  points  are 
connected  into  a  tree  via  the  partial  matches.  The  full  matches  (permutations)  form  the  leaves  of 
the  tree.  The  hope  is  that,  although  the  tree  contains  k!  leaves,  most  of  the  branches  can  be  pruned 
(removed  from  further  consideration)  early.  The  elegant  approach  of  Baird  also  employs  the  tree- 
pruning  technique  and  claims  to  be  computationally  superior  to  the  other  methods.  The  inherent 
problems  associated  with  the  tree-pruning  methods  are  explained  in  [13]. 

Image  registration  is  often  formulated  as  a  “least-squares”  problem.  Partial  results  toward  the 
solution  of  the  general  problem  have  been  discovered  and  rediscovered  several  times.  Perhaps  the 
oldest  paper  on  the  subject  is  green  [8],  1952.  Long  before  various  computer  vision  problems  became 
the  most  pressing  open  problems  of  the  robotics- computer  science  of  today,  B.  Green  gave  a  solution 
to  the  optimal  rotation  subproblem.  Green’s  result  was  later  improved  by  shonemann  [11].  Both 
papers  appeared  in  Psychometrics  as  the  research  focused  on  factor  analysis  rather  then  on  image 
registration.  Clearly  unaware  of  the  classica'  results,  faugeras  and  hebert  [5]  (in  a  more  general 
work)  gave  solutions  to  the  two-  and  three-dimensional  rotation  subproblems  by  a  method  different 
from  that  of  Green  and  Shonemann.  Still  later,  arun,  HUang,  and  blostein  [3]  rediscovered  the 
original  method  in  the  image  registration  setting.  The  latter  two  papers  also  give  partial  results 
on  the  optimal  translation  subproblem.  All  these  results  implicitly  or  explicitly  assume  that  the 
one-to-one  matching  of  points  is  either  fixed  or  known.  These  amd  other  aspects  of  the  general 
least-squares  problem  are  treated  in  21KAN  and  SILBERBERG  [13].  An  O(ib^)  approximation  solution 
method  to  the  general  problem  is  also  there  given. 

For  other  formulations  of  the  image  registration  problem,  see  BAIRD  (4),  ALT,  mehlhorn,  wagener, 
and  WELZL  [l],  and  ARKiN,  MITCHELL,  and  ZIKAN  [2].  A  Strongly  polynomial  algorithm  which  solves 
the  Baird  s  formulation  can  be  (essentially)  found  in  [1].  In  [2]  a  strongly  polynomial  algorithm  to  a 
more  general  problem  is  given.  Most  natural  formulations  of  the  image  registration  problem  can  be 
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at  loast  approximately  cast  within  this  framework;  the  “least-squares”  and  Baird  formulations  are 
no  exceptions.  If  we  can  extrapolate  from  the  experience  of  other  branches  of  applied  mathematics, 
however,  then  we  can  see  that  a  good  specialized  algorithm  for  the  “least-squares”  registration 
formulation  is  likely  to  be  computationally  superior  to  other  serious  image  registiation  approaches. 
We  utilize  parametric  linear  programming  to  provide  such  an  algorithm  here. 

The  following  notions  (and  the  associated  notatbn)  from  complex  plane  geometry,  mathematical  pro¬ 
gramming.  and  metric  space  theory  provide  a  convenient  formal  framework  for  the  two-dimensional 
image  registration  problem. 

is  a  two-dimensional  vector,  then  let  *  =  xi  -P  ixj  =  [||x||,0c]  be  the  corresponding 
complex  number,  where  ||x||  =  ||x||j  is  the  magnitude  (two-norm).  Or  is  the  argument,  and  [||x||,^x] 
is  the  polar  representation.  Let  JT  denote  the  complex  r.oryugate  of  x  and  recall  that  in  polar 
coordinates  xy  =  (||x||  i|y||.6r  +  xy  =  (||x||  •  ||y||,9x  -  6»)]  The  scalar  product  of  x  and  y 

is  x  o  y  =  ri.i/i  -(-  xjyi  -  ||x!|||y|| cos(0x  -  tfy),  the  “cross”  product  is  x  x  y  =  ||x||||y|| sin(tfx  -  0^). 
Recall  that  xy  =  x  o  y  -p  ix  x  y  and  that  =  x  o  y.  The  scalar  product  is  also  called  the  “inner” 
or  the  “dot”  product.  To  prevent  possible  ambiguity  in  mathematical  formulas,  let  the  notation 
(•,  )  denote  the  inner  product.  Therefore  (x,y)  =  xc  y  for  complex  numbers,  (x,y)  =  x’'j/  for  real 
^  dimensional  vectors,  and  if  A  and  B  Me  real  n  x  it  matrices,  then  let  {A,B)  =  trace(wdiJ^). 

Assume  that  A  and  D  are  k-diiiKnsional  vectors  of  complex  numbers  and  define 

k 

(A,B)  =  ^.jE3,  (l-l) 


where  aj  and  bj  are  the  j-th  components  of  the  vectors  A  and  D,  respectively.  This  is  the  complex 
inner  product. 

Consider  the  complex  number  e**  =  coe(^ )  -p  i  8in(tf)  Counterclockwise  rotation  of  two-dimensional 
vectors  by  the  angle  0  can  be  represented  by  multiplication  of  the  corresponding  complex  numbers 
by  e'*.  If  r  =  [sr.^r],  then  the  multiplication  by  r  corresponds  to  the  counterclockwise  rotation  by 
Or  and  scaling  by  Sr  >  0.  Let  a,  denote  the  argument  of  the  complex  number  at  and  /?;  be  the 
argUFTK-nt  of  bp  Define 

d.,(r)  =  ||a,  -  rbjll*  =  ||a,||^  -P  llrl^lbjH^  -  2{a,.rbj) 

(1  -2) 

“  ll«ill^  +  -  2«r||ai||||bj!|c06(rt,  -0J  -  Or) 


for  all  distinct  pairs. 


I 


3 


Consider  llie  parainelric  linear  progranrs 


min  min  >>  d,,(r)  *,, 
r€C'>  X 

1=1  j  =  i 
t 

subject  to  ^  x,j  =  1,  for  all  j  =  1,2, . . . ,  t 

■=i 

k 

^  r,y  =  1,  for  all  i  =  1, 2, . . . ,  t 


f,,  >  0,  for  all  I,  j  =  1, 2, . . . ,  t, 


(1-3) 


mill  min>  >  ■  Xj, 

!|r|l  =  l 

1=1  ;=1 
k 

ubject  to  ^^Xij  =  1,  for  all  ji  =  1,2 . k 

i 

k 

^  z,j  =  1,  for  all  i  =  1 , 2, . . . ,  t 
j  =  i 

x,^  >  0,  for  all  i,  j  =  1, 2, . . . ,  k. 


(1-4) 


These  problems  arise  naluraJly  in  the  context  of  “least-squares”  image  registration,  zikan  »nd 
siLDERBERG  [13].  If  (r*,x*)  is  an  optimal  solution  of  (1-3),  then  point-to-point  matching  can  be 
recovered  from  r*,  and  rotation  and  scaling  from  r' 

It  is  natural  to  view  an  n-dimensional  image  consisting  of  k  points  as  a  t  x  n  matrix  with  rows 
corresponding  to  the  appropriate  image  points.  If  the  image  is  two-dimensional,  then  it  is  also 
convenient  to  associate  with  the  image  a  Jfc-dimensional  complex  vector  A  where  each  coordinate  ai 
corresponds  to  one  image  point.  In  this  paper  we  choose  the  latter  formalism.  The  association  of 
a  vector  (mv'rix)  to  an  image  is  not  unique;  it  depends  on  the  ordering  of  the  points.  In  general, 
consider  the  complex  t-dimensional  vector  space  C* .  Assume  that  r  and  t  are  complex  numbers.  The 
image  as.sociated  with  the  vector  A  6  C*  can  be  translated  by  A  0  t  =  {ai  -I-  t,aj  -b  t, . .  . ,  -|-  t) 

and  rotated  and  seeded  by  rA.  AfTinc  transformation  of  image  A  can  then  be  written  as: 


(r,t){A)  =  r(A0t). 


(1-5) 


A  relabeling  of  image  points  by  permutation  ir  acts  on  the  image  by  )r(A)  =  where  Pw  is  the 

permutation  matrix  corresponding  to  x.  Two  images  are  to  be  considered  as  being  eqniva/c/it  if  and 
only  if  they  differ  by  a  registration  and  relabeling  only. 


A  =  (r,rr,t)D  =  r/\(D  0  t). 


(1-6) 
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Recall  that  ||A||  =  (A,A)^  is  the  euclidean  (FVobenius)  norm  of  the  vector  A.f  Assume  that 
vector  A  corresponds  to  an  image  which  has  been  translated  rotated  and  scaled.  Assume  also  that 
noise  has  been  added  to  each  point.  Denote  the  original  image  by  B.  If  we  predict  the  unknown 
transformation  from  the  optimal  solution  of  the  problem 

min  ||A  -  (r,»,t)B||,  (1-7) 

then  we  have  the  “least-squares”  estimate,  in  spirit  analogous  to  least-squares  estimates  in  other 
brtnehes  of  science.  The  matching  and  transformation  give  the  least  sum  of  the  squared  distances 
between  the  matched  points.  In  [13]  we  show  that  (1-7)  essentially  reduces  to  (1-3)  or  (1-4).  Thus 
these  parametric  programs  are  important  from  the  practical  standpoint. 

2.  On  Frobenius  Norm  Image  Registration 

Many  practical  aspects  of  using  the  Frobenius  norm  based  criteria  for  n-dimensional  image  regis¬ 
tration  were  addressed  in  zikan  «nd  siLBERBi::ttG  [13].  Most  importantly,  the  issue  of  missing  and 
spurious  points  was  discussed.!  In  this  paper  we  restrict  our  discussion  to  the  case  where  no  missing 
or  spurious  points  occur.  The  extension  of  our  results  to  the  unequal  cardinality  case  would  be  done 
analogously  to  the  development  in  |13]. 

On  the  theoretical  side,  it  is  shown  in  [13]  that  the  optimal  translation,  independent  of  rotation, 
permutation  or  scaling,  always  translates  the  center  of  mass  (centroid)  of  B  to  the  center  of  mass  cf 
A.  Let  A  and  B  be  n  x  I:  real  or  complex  matrices,  and  a  and  0  be  the  respective  centers  of  mass 
(centroids)  of  their  row  vectors.  (For  us,  n  =  I  and  the  coordinates  are  viewed  as  image  points.) 

Theorem  1.  If  A  and  B  are  n  x  k  matrices  owr  the  real  or  complex  Reid  V ,  then  a  —  0  solves  the 
translation  problem 

mmJIA  -  (B®t)||r.  (2-1) 

Proof:  A  simple  proof  from  first  principles  can  be  found  in  (13,  Section  3.1).  o 

The  centroids  of  the  images  can  be  superimposed  and  conveniently  identified  with  the  origin  of  the 
cartesian  coordinate  system. 

It  is  also  argued  in  [13]  that  the  scaling  parameter  can  be,  and  perhaps  should  be,  estimated  before 
one  begins  the  computation  of  the  optimal  rotation  and  matching.  This  result  is  enhanced  later  in 

t  In  the  case  of  n-dimensional  images,  where  the  k  x  n  real  matrices  take  the  place  of  vector  A, 
we  use  the  Frobenius  matrix  norm  ||A/||j.  =  (trace(Af M'^))^  =  (jV/,Af)4.  The  Frobenius  metric  is 
the  same  as  the  euclidean  norm  when  M  is  viewed  as  a  kn-dimensiond  vector.  If  Af  is  a  complex 
matrix,  then  HAfHj-  —  tracefAf  Af*),  where  Af*  =^7^. 
t  The  presence  of  (missing  and)  spurious  points  leads  to  the  problem  of  registering  images  of  different 
cardinalities. 
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tliis  paper,  when  it  is  shown  that  positive  scaling  has  no  effect  on  optirnal  rotation  and  matching. 
If  the  matching  is  fixed,  then  the  optimal  rotation  and  scaling  of  a  two-dimensional  registration 
problem  can  be  recovered  with  the  help  of  the  following  theorem. 


Tlicorein  2.  Let  A  and  D  be  k-JiniensioriAl  complex  vectors,  such  that  (A,D)  ^  0.  If 

it  (A,B) 

■  ll(A,B)|| 

and 

^  _  ll(A.D)|| 

IIDIP  ■ 

then  o'*  solves  the  problem  min«i>  1|A  —  c**D||.  and 

.s  _  (A.D) 

-  ||D|p 

solves  miiir  1|A-  rD||. 

Proof:  For  a  proof  of  the  theorem  consult  (13,  Section  4,1  and  Appendix  1],  □ 


(2-2) 


(2  -  3) 


(2-4) 


.\nalogous  results  for  rotation  and  scaling  in  higher  dimensions  can  be  obtained,  [13,  Section  4].  The 
three-dimensional  problem  can  be  solved  with  the  help  of  the  “algebra”  of  quaternions,  faugeras  »nd 
HEBERT  [5],  and  the  general  n-dimensional  problem  with  the  help  of  the  singular  value  decomposi¬ 
tion.  GREEN  [8],  and  smonemann  [11],  The  optimal  rotation  problem  is  known  in  literature  as  the 
orthogonal  Procrustes  problem,  COEUB  •nd  van  loan  (7),  after  the  villainous  son  of  Poseidon.* 

The  case  that  is  not  explicitly  covered  by  the  theorem,  i.e.  (A,D)  =  0,  may  be  treated  by  one  of 
the  standard  lexicographical  methods,  zikan  and  sieberberg  (13). 

In  fl3]  one  finds  an  approximation  scheme,  based  on  parametric  linear  programming,  which  solves 
(1-4)  (and  consequently  (1-3))  within  aspecified  erior  in  O(t^)  worst-case  computational  complexity. 
The  availability  of  “hot  start"  at  each  step  of  the  parametric  “sweep"  luakes  the  approach  practically 
attractive.  In  this  paper  we  enhance  the  parametric  algorithm.  The  new  algorithm  enables  us  to 
solve  (1-3)  and  (1-4)  exactly  and  (in  fact)  faster  than  by  the  original  approximation  method.  In  the 
process  of  developing  the  algorithm  we  also  enrich  the  theory  of  the  least-squares  formulation  of  the 
image  registration  problem. 


*  In  Greek  mythology,  Prociustes  forced  travelers  to  fit  into  his  bed  by  stretching  their  bodies  or 
cutting  off  their  legs. 
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3.  Bilineitr  Functionals  and  the  Parametric  Linear  Programming  Problems 


Define 

A  =  |o  <  r  €  ~  ***  ”  1.2,...  ,*|  .  (3-1) 

Recall  the  classical  result  due  to  Birkhoff  (1946)  that  the  vertices  of  X  correspond  to  permutations 
of  k  elements.  The  vector  r  is  a  vertex  of  -V  if  and  only  if  there  exists  a  permutation  a  such  that 
x,j  =  1  whenever  t  =  ir(j),  and  Zi^  =  0  otherwise.  The  bases  of  X  correspond  to  the  spanning  trees 
of  the  (complete)  bipartite  graph  Bkk 

Recall  the  definition  (1-2)  of  d,,(r).  Assume  a  lexicographic  ordering  of  the  {i,j)  pairs  and  let 


d(r)  =  (dij(r)j  denote  the  resulting  real  b'-dimensional  vector, 
be  written  as 

Note,  that  (1-3)  and  (1-4)  can  now 

min  min(d(r),x), 

r  A' 

and 

(3-2) 

trjin  min(d(r).*). 

l|r||=>*€A 

Define  the  natural  functionals, 

{3-o> 

F(r,x)  =  (d(r).i) 

and 

(3-4) 

F(r)  =  min  F(r,x), 

'  '  *€A  ' 

(3-5) 

asscKiated  with  (3-2)  and  (3-3).  In  addition  to  X,  it  is  convenient  to  introduce  analogues  of  d(r), 
F(r,x),  and  F(r).  Thus,  let  us  define  the  auxiliary  functions 

^.>(«')  =  -(“i.'bj), 

(3-6) 

G(r,x)  =  (6(r),x), 

and 

(3-7) 

G(r)  =  ininG(r,x). 
*€X 

(3-8) 

Note  that  C{  .  )  is  a  biiinear  funct/onaJ.  Also  note  that 


<f.,(r)=  ll‘‘<ll*  +  IMl’ll»>jll’  +  2«.;(r).  (3-9) 

and 

F(r,  X)  =  l|A|p  -k  ||r|p||B||^  -k  2C(r,x).  (3  -  10) 

Lemma  1.  Junctions  F(r)  and  G(r)  are  related  by 

F(r)  =  II AH'  -k  ||r||'||B||*  +  2G(r).  (3-11) 
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Proof:  Pqiialioii  (3-11)  results  from  these  idenlities: 


F(r)  =  ii»in  F(r,x)  =  min(<f(r),x) 

xCA’ 

k 

~  llrll^ll^ijll*  +  2<,j(r)}  x,j 

*.>=M 

^  i=i;=i  '  ;=ii=i  ij=i 

=  ->■  2min{i(r).r) 

>=1  J=1  >=1  •=! 

=  l|All^-H|r|niB||’-h2G(r), 

where  the  first  three  identities  follow  from  (3-5),  (3-4),  and  (3-9)  respectively,  and  the  last  relation 
holds  since  x,j  =  52*_i  i,j  =  I,  (3-1).  Simple  algebraic  manipulations  complete  the  proof  o 


Thooreiii  2.  The  complex  number  «e'**  solves  niin||r||=i  C(r)  if  and  only  if  it  solves  min||,||_i  F{r). 
Proof:  If  the  norm  of  r  is  fixed,  then  by  Lemma  1,  P(r)  and  2G(r)  differ  by  a  constant  only,  o 


Above  all,  we  are  interested  in  the  case  s  =  1.  Theorem  2  implies  that  if  the  scaling  factor  is 
known,  then  the  least-squares  problem  associated  with  the  penalty  function  d  is  equivalent  to  the 
■maximal-projection”  problem  associated  with  the  penalty  function  6. 

Lonima  3.  Assume  that  s  >  0  is  a  nonnegative  sca/v.  If  r*  €  -V  solves  the  linesr  progrvn  (3-8) 
associated  with  G(r),  then  it  solves  the  linesr  program  associated  with  G(sr).  The  optimal  solutions 
are  related  by 

G(sr)  =  sG(r).  (3-13) 

Proof;  Since  G(r,  i)  is  a  bilinear  functional, 

(«(sr),x)  =  {s6(r),r)  =  s(«(r),r). 

If  s  =  0,  then  the  result  is  immediate.  If  s  >  0,  then  the  lemma  readily  follows  because  positive 
scaling  of  the  objective  function  of  a  linear  program  does  not  change  its  optimal  solution  set.  o 
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TliiHireiii  4.  Assuiiio  t/i«C  x’  €  .V  is  a  vertex  so/ution  to 


nun  minGtr,  x) 

||r||Etr€X 


(or  sonie  positive  s  Assume  that  x*  is  the  permutation  aasoriated  with  z*,  and  Jet 


o‘**  = 


Then  {o**’,x*}  selves  the  prahicint 


(a)  min  minGJr.z), 
|M|=lrCX 


and 

(b)  inin  niinFtr.z), 

11,11=1  »€X 

while  {r*,x*}  solves  the  problem 

(c)  min  inin  Ffr.x) 

T  rfX 


(3-14) 

(3-15) 

t3-lC) 

(3-17) 

(3-18) 


Proof:  Wc  assume  that  there  exists  c'*  5U<-Ii  that 

G(o‘»'.x*)  <  G(c**,x)  (3-19) 

for  all  c'*  and  x,  Conrcquently,  by  Theorem  3.2, 

(■'{e**',w’)  <  F{e'*,ir)  (3-20) 

for  all  e'*  and  x.  But  we  know  from  Theorem  2.2  that 

F(e‘*’,x-)  <  F(e‘',x*).  (3-21) 

for  all  c**.  It  follows  from  (3-20)  and  (3-21)  that 

^(c‘*',x*)  <  f(c'*',x-)  <  F{e‘*,x)  (3-22) 

for  all  e**  and  x.  This  establishes  that  {o'*  ,x*)  solves  the  problem  (b).  Invoke  the  Theorem  3.2 
(again)  and  from  (3-22)  obtain 

<  G(e“.ir)  (3-23) 
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f('r  all  *>'*  Aiid  IT,  which  pstablishcs  that  {c'**,/*)  sulvet  the  problem  (a)  Finally,  for  all  nonncgativc 

.V  >  0  wc  have 

=  «G(o‘**.»*)  <  aC{c'*,ir)  =  0’(«o‘*,ir)  (3-24) 

by  I.eiiiiiia  3  3  Invoke  Theorems  2.2  and  3.2  together  and  obtain 

<  F(sc'*‘,w')<  (3-25) 


Since  s  is  arbitrary,  we  have  completed  the  proof  that  {r*,i*)  solves  the  problem  (c)  o 


l.et  /east -$(jn ares  regisirarion  problems  be  a  generic  name  for  ti>e  three  problem^  of  (3- 18)  We 
ii:\v<  established  a  strong  lelation  between  the  functionals  G(  ,  )  and  F(-,  •).  The  bilinear  slrnclnre 
of  f>'(  ,  )  provides  us  witn  anotlier  special  property  of  the  optimal  solution  sets.  After  we  establish 
tins  last  nontrivial  theoretical  property,  we  can  give  the  promised  algoriclim.  For  a  fixed  0,  define 
the  corresponding  complex  plane  equivalence  relation  by 


r(»  if  and  only  if 

G(r''',r-‘)  <  G(r'*’.r)  for  all  j  6  A'  <=>  G(r^*\x*)  <  G(r**\x)  for  all  x  6  A. 


(3  -  26) 


Tlu'orcin  5.  ReUtior)  induces  a  convex  cuniceJ  subdivision  of  the  complex  plane. 

Proof:  Assume,  that  (3  26)  holds  for  some  pair  r**’  and  Then  by  Lemma  3  3 

G(5r“>,  J*)  =  ,'!G(r'*>.T*)  <  *G(r‘*>.x)  =  G(.sr“>,x)  (3  -  27) 

for  all  s  >  0,  and 

G(Ar‘*\T’)  +G((1  -  A)r'’>,i*)  = 

AG(r'*',x*)  +  (1  -  A)C;(r<*»,x*)  <  AG(r'*>.x)  +  (1  -  A)G(r‘*>.x)  (3  -  28) 

=  G(Ar“».x)  +  G((l-A)r<»,x), 

for  all  0  <  A  <  1.  0 


'1  lie  subdivision  induced  by  F  is,  of  course,  identical  to  the  one  induced  by  G  and  consists  of  a  finite 
numb,  r  of  cells.  Thus  the  subdivision  is  the  union  of  closures  of  its  ''two-diiiiensionBl'’  colls.  (If  we 
\  u  w  {’  as  a  two-dimensional  real  plane.)  For  each  two-dimensional  cell,  <r,  we  arbitrarily  choose  a 
represenratne  xftr)  from  the  vertices  of  the  associated  solution  set  of  niiii,.(A'  G(r„,x),  where  r„  is 
any  element  of  <t.  Let  n(<r)  be  the  pcrmut,ation  corresponding  to  xfrr) 
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Proposition  6.  The  set  of  representatives  contains  at  least  one  optimal  matching  associated  with 
the  least-squares  registration  problems  (3-18). 

Proof:  The  result  follows  from  the  fact  that  the  relation  G(r,r*)  <  G{r,  x)  for  all  r  €  extends  by 
continuity  to  the  boundary  of  <t.  □ 

Figure  3.1  exhibits  a  finite  coniceil  subdivision  of  the  plane. 


r 

2 


r 

1 


Figure  3.1  .4  conical  subdivision  of  C  and  a  piecewise-linear  path  closed  around  origin 

Note,  that  any  closed  curve  around  the  origin  intersects  all  two-dimensional  regions  of  the  subdivi¬ 
sion.  In  particular,  this  is  true  of  all  piecewise-linear  curves  around  the  origin. 
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4.  The  Algorithm 


The  algorithm  to  solve  the  least-squares  registration  problems  (3-18)  has  four  bsisic  steps: 

ALGORITHM 

STEP  1:  Choose  a  closed  ptecewtse-linear  curve  in  the  plane  fyj,  which  contains  origin  in  its 
interior. 

STEP  2:  Choose  a  starting  point  r°  €  7  and  solve  the  assignment  problem  associated  with  G(r°). 

STEP  3:  Parametrically  compute  the  optimal  solutions  of  the  assignment  problems  associated  with 
G(y(t))  and  collect  all  distinct  locally  optimal  permutations. 

STEP  4:  Use  Theorem  2.2  to  find  the  optimal  rotation  and  scaling  corresponding  to  each  permu¬ 
tation  obtained  in  Step  S.  The  best  overall  solution  solves  the  registration  problems  (3-18). 

Let  us  briefly  remark  on  each  individual  step. 

Step  J  The  unit  square  of  Figure  3,1  is  a  convenient  choice  for  7.  Note,  that 

=  (4-1) 

and  that 

^,(>)  =  aixEj,  (4-2) 

the  “dot”  and  “cross”  products  of  ai  and  Ej  respectively.  If  r  =  rj  ir^,  then 

iij{r)--riS,j{l)  +  rj6,,{i).  (4-3) 


Consequently,  the  costs  on  the  unit  square  of  Figure  31  are  easy  to  construct. 

Step  2.  It  is  convenient  to  choose  a  cornerpoint  of  the  unit  square  as  a  starting  point,  for  instance, 
r®  =  1  -I-  i.  The  Hungarian  method  can  be  used  to  solve  G(1  -I-  i)  in  0(1^)  worst-case  “time”. 

Step  3  The  parametric  version  of  the  simplex  algorithm  specialized  to  the  transportation  problem 
can  be  used  to  perform  this  step.  Since  the  vertices  of  X  (3-1)  are  degenerate,  few  blocked  pivots 
may  be  performed  between  successive  permutations.  It  is  believed  that  the  overall  number  of  pivots 
is  quadratic  in  k,  however,  this  question  is  still  open.f 

t  Consider  the  conical  subdivision  of  generated  by  the  nonnegativity  constraints  in  the  circu¬ 

lations  subspace,  papadimitriou  and  STEIGLITZ  (10),  KENNINGTON  and  HELGASON  (^9).  The  Complexity 
of  our  algorithm  directly  depends  on  the  number  of  cones  intersected  by  a  line  through  this  conical 
subdivision. 
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Step  4.  In  practice,  Steps  3  and  4  are  merged  so  that  not  all  relevant  matchings  need  be  stored 
in  memory.  Note  (Theorem  2.2)  that  for  each  permutation  we  mainly  need  to  compute  (A,ir(B)), 
(l-l).  This  requires  0{k)  arithmetic  operations. 
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Imife  regiAtratioD  involves  estimating  how  one  set  of  n-dimensional  points  ia  rotated,  scaled,  and 
translated  into  a  secood  set  of  n-dimensional  points.  1q  practice,  n  is  usually  2  or  3.  We  give  an  exact 
algorithm  to  solve  the  "least-squares"  formulation  of  the  two-dimensional  registration  problem.  The 
algorithm,  which  ia  based  on  parametric  linear  programming,  can  be  viewed  as  a  refinement  of  the 
0(k*)  approximation  method  proposed  by  Zikan  and  Silberberg  (13).  The  approach  can  be  extended 
to  handle  registration  of  images  of  different  cardinalities. 
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